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With particular attention to the recently postulated introduction of a non-extensive generalization 
of Boltzmann-Gibbs statistics, we study the long-term stellar dynamical evolution of self-gravitating 
systems on timescales much longer than the two-body relaxation time. In a self-gravitating A''-body 
system confined in an adiabatic wall, we show that the quasi-equilibrium sequence arising from the 
Tsallis entropy, so-called stellar polytropes, plays an important role in characterizing the transient 
states away from the Boltzmann-Gibbs equilibrium state. 



Introduction. — The long-term evolution of self- 
gravitating many-body system is an old problem with 
a rich history in astrophysics. The problem, in nature, 
involves the long-range nature of attractive gravity and 
is fundamentally connected with statistical mechanics 
and thermodynamics. Historically, the important conse- 
quence from the thcrmodynamical arguments has arisen 
in 1960s, known as the gravothermal catastrophe, i.e., 
thermodynamic instability due to the negative specific 
heatQ. Originally, the gravothermal catastrophe has 
been investigated in a very idealized situation, i.e., a stel- 
lar system confined in a spherical cavity 0- Owing to the 
maximum entropy principle, the existence of an unstable 
thermal state has been found from the standard analysis 
of statistical mechanics with a particular attention to the 
Boltzmann-Gibbs (BG) entropy: 



j d^xd^v f{x, v) In f{x, v), 



(1) 



where f(x,v) denotes the one-particle distribution func- 
tion defined in phase-space {x,v). 

Since 1960s, the standard approach using the BG en- 
tropy has dramatically improved our view of the late- 
time phase of the globular cluster as a real astronomical 
systemjSjj, however, from a thermodynamic point of view, 
peculiarities of the thermal property in self-gravitating 
systems such as negative specific heat, as well as the non- 
equilibrium properties away from the BG state have not 
yet been understood completely. 

In this Letter, aiming at a better understanding of the 
(non-equilibrium) thermodynamic properties of stellar 
self-gravitating systems, we present a set of long-term iV- 
body simulations, the timescale of which is much longer 
than the relaxation time. With a particular emphasis 
to the recent application of the non-extensive generaliza- 
tion of BG statistics, we focus on the stellar dynamical 
evolution in an isolated star cluster before self-similar 
core-collapse*^. We show that the quasi-equilibrium se- 
quence arising from the Tsallis entropy 5j plays an impor- 
tant role in characterizing the non-equilibrium evolution 
of a self-gravitating system. 



N-body simulations. — The A'^-body experiment con- 
sidered here is the same situation as investigated in clas- 
sic papers (12, see also Ref.Q). That is, we confine the 
N particles interacted via Newton gravity in a spherical 
adiabatic wall, which reverses the radial components of 
the velocity if the particle reaches the wall. Without loss 
of generality, we set the units a,s G — M — r^ — 1, where 
G is gravitational constant, M and is the total mass 
of the system and the radius of the adiabatic wall, re- 
spectively. Note that the typical timescales appearing in 
this system are the free-fall time, Tg = (Gp)~^/^, and 
the global relaxation time driven by the two-body en- 
counter, Ti-eiax = (O.liV/ In A^)Tfi 0, which are basically 
scaled as Tff ^ 1 and Trciax ~ O.liV/lniV in our units. 
To perform an expensive N-body calculation, we used a 
special-purpose hardware, GRAPE-6, which is especially 
designed to accelerate the gravitational force calculations 
for coUisional iV-body systems With this implemen- 
tation, the 4th-order Hermite integrator with individual 
time step 8j can work efficiently, which is suited for prob- 
ing the relaxation process in denser core regions with an 
appropriate accura cy. We adopt the Plummer softened 
potential, 4> = l/y/r^ + e^ with a softening length e of 
1/512 and 1/2048. 

In the present numerical simulation, the choice of the 
initial condition is an important step. Here, we set the 
initial distribution to the stationary state of Poisson- 
Vlassov equation, i.e., dynamical equilibrium for a spher- 
ical system with isotropic velocity distribution. Accord- 
ing to the Jeans theorem JJ, the one-particle distribution 
function f{x, v) can be expressed as a function of specific 
energy, e = v'^ /2-\-^{r) with r and $ being the radius and 
the gravitational potential. Then keeping the energy and 
the mass constant, the thermal equilibrium of ordinary 
extensive statistics derived from the maximum entropy 
principle of the BG entropy reduces to the exponen- 
tial distribution, so-called isothermal distribution given 
by /(e) oc e~^^, which effectively satisfies the equation 
of state, P(r)c(. p(r), where P{r) is pressure and p{r) is 
mass density 111, • 

On the other hand, as another possibility, one considers 
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TABLE I: Initial distributions and their evolutionary states 



run # initial distribution parameters # of particles transient state final state 

A stellar polytrope(n = 3) D = 10, 000 2,048 stellar polytrope collapse 

Bl stellar polytrope(n = 6) _D = 110 2,048 stellar polytrope collapse 

B2 stellar polytrope(n = 6) D = 10, 000 2,048 stellar polytrope collapse 

CI Hernquist model a/re = 0.5 8,192 stellar polytrope collapse 

C2 Hernquist model a/re = 1-0 8,192 none isothermal 




10=^ 10= 10" 10= 

D(=Pc/Pa) 

FIG. 1: Equilibrium sequences of stellar polytrope and 
isothermal distribution (n = oo) in the energy-density con- 
trast relation, A = -TeEjiGM ) v.s. D = pc/pe- The thick 
arrows denote the evolutionary tracks in each simulation run. 



the extremum state of Tsallis' non-extensive entropy 0: 



S,^- I d^xd^v [{fix, v)r - fix, v)]/il - q), (2) 



which might be of particular importance in describing 
the quasi-cquilibrium state away from the BG stateQ. 
In this case, the maximum entropy principle leads to the 
power-law distribution, /(e) oc ($o - e)^/'^''"^^ [lo l ITU 
lillil, referred to as the stellar polytrope^ . It satisfies 
the polytropic equation of state, P(r) oc /9(r)^"'"^/", and 
the polytrope index n is related to the g-parameter as 
n = — 1) -1-3/20. Provided the polytrope index n, 
the equilibrium structure can be determined by solving 
the Lane-Emden equation and using this solution, 
the relationship between the dimensionless energy A = 
—reE/iGM'^) and the density contrast D = Pc/pe, the 
core densit y d ivided by the edge density, can be drawn 
(see Fig^ Note that the limit n ^ oo (or g ^ 1) 

corresponds to the isothermal distribution derived from 
the BG entropy iQ}. 

Table ^summarizes the list of the five simulation runs. 
A more systematic study of the systems with several ini- 
tial conditions is now in progress and the details of the 
results will be reported elsewhere. In Table in addi- 



tion to the stellar polytropic initial state, we also consider 
the Hernquist model(l7l|. which was originally introduced 
to account for the empirical law of observed elliptical 
galaxies Q. 

Results. — It has been recently discussed in pTllT^IT^ 
that the thermodynamic structure of the stellar poly- 
tropic distribution can be consistently characterized by 
the non-extensive framework of the thermostatistics. Ac- 
cording to their results, the stellar polytrope confined in 
an adiabatic wall is shown to be thermodynamically sta- 
ble when the polytrope index n < 5. In other words, 
if n > 5, a stable equilibrium state ceases to exist for 
a sufficiently large density contrast D > Z^crit, where 
the critical value Dcrit given by a function of n is de- 
termined from the second variation of entropy around 
the extremum state of Tsallis entropy, 6^Sq = O jlll [l3l |. 
The dotted line in Fig represents the critical value Dcrit 
for each polytrope index, which indicates that the stellar 
polytrope at low density contrast D < Dcrit is expected 
to remain stable. Apart from the BG state, one might 
expect that the stellar polytrope acts as a thermal equi- 
librium. 

Of course, this naive expectation is not correct at all. 
Indeed, the numerical simulations reveal that the stellar 
polytropic distribution gradually changes with time, on 
the timescale of two-body relaxation. Further, it seems 
that the gravothermal instability appears at relatively 
smaller values of D than the predicted one, -Dcrit, which 
might be partially ascribed to the non-stationarity of 
the background stellar polytrope. Physically, the core- 
collapse phenomenon due to the gravothermal catas- 
trophe follows from the decoupling of the relaxation 
timescales between the central and the outer parts, whose 
behavior sensitively depends on the physical property of 
heat transport [isj. In a rigorous sense, the thermody- 
namic prediction might lose the physical relevance, how- 
ever, focusing on the evolutionary sequence, we found 
that the transient state starting from the initial stellar 
polytrope can be remarkably characterized by a sequence 
of stellar polytropes (run A, Bl and B2). This is even 
true in the case starting from the Hernquist model (run 
CI). 

Let us show the representative results taken from run A 
(Fig|2J). Fig[2Ia) plots the snapshots of the density profile 
p(r), while Figl^Jb) represents the distribution function 
fie) as a function of the specific energy e. Note that 
just for illustrative purpose, each output result is arti- 
ficially shifting to the two-digits below. Only the final 



FIG. 2; Results from simulation run A. (a) Snapshots of density profile p{r). (b) Snapshots of one-particle distribution function 
/(e). (c) The time evolution of the polytrope index for run A with and without the boundary of adiabatic wall. 





FIG. 3: Evolution of one-particle distribution function in other models, (a) Run Bl. (b) Run CI. (c) Run A removing the 
adiabatic wall. 



output with T ~ 400 represents the correct scales. In 
each figure, solid lines mean the initial stellar polytrope 
with n = 3 and the other lines indicate the fitting results 
to the stellar polytrope by varying the polytrope index 
n[l9|. Note that the number of fitting parameters just 
reduces to one, i.e., the polytrope index, since the total 
energy is well-conserved in the present situation. Fig[21 
shows that while the system gradually deviates from the 
initial polytropic state, the transient state still follows 
a sequence of stellar polytropes. The fitting results are 
remarkably good until the time exceeds T ~ 400, corre- 
sponding to ISTiciax- Afterwards, the system enters the 
gravothermally unstable regime and finally undergoes the 
core-collapse. 

Now, focus on the evolutionary track in each simula- 
tion run summarized in the energy-density contrast plane 
(Fig^ , where the filled circle represents the initial stellar 
polytrope. Interestingly, the density contrast of the tran- 
sient state in run A initially decreases, but it eventually 
turns to increase. The turning-point roughly corresponds 
to the stellar polytrope with index n ~ 5 — 6. Note, how- 
ever, that the time evolution of polytrope index itself is 
a monotonically increasing function of time as shown in 
FiglUc), apart from the tiny fluctuations. This is indeed 
true for the other cases, indicating the Boltzmann H- 
theorem that any of the self-gravitating systems tends to 
approach the BG state. A typical example is the run C2, 
which finally reaches the stable BG state. However, as 



already shown in run A, all the systems cannot reach the 
BG state. Fignindicates that no BG state is possible for 
a fixed value A > 0.3350, which can be derived from the 
peak value of the trajectory. Further, stable stellar poly- 
tropes cease to exist at high density contrast D > Z^crit- 
In fact, our simulations starting from the stellar poly- 
tropes finally underwent core-collapse (runs A, Bl and 
B2). Though it might not be rigorously correct, the pre- 
dicted value I?crit provides a crude approximation to the 
boundary between the stability and the instability. 

Fig|3| plots the snapshots of the distribution function 
taken from the other runs. The initial density contrast 
in run Bl (Figl^Ja)) is relatively low {D = 110) and 
thereby the system slowly evolves following a sequence 
of steflar polytropes. After T = 2000 ~ 74Trciax, the 
system begins to deviate from the stable equilibrium se- 
quence, leading to the core-collapse. Another noticeable 
case is the run CI (Fig|3Ib)). The Hernquist model 
as initial distribution of run C has cuspy density pro- 
file, p{r) cx l/r/(r + a)^, which behaves as p cx at 
the inner part 17]. The resultant distribution function 
/(e) shows a singular behavior at the negative energy 
region, which cannot be described by the power-law dis- 
tribution. Soon after a while, however, the gravothermal 
expansion^ takes place and the flatter core is eventually 
formed. Then the system settles into a sequence of stel- 
lar polytropes and can be approximately described by the 
stellar polytrope with index n = 20 for a long time. Thus, 
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as long as the stellar system confined in an adiabatic wall 
is concerned, the stellar polytrope can be regarded as a 
quasi-attractor and a quasi-equilibrium state. 

Of course, these remarkable featm'cs could be an out- 
come in a very idealized situation and one suspects that 
quasi-equilibrium state of stellar polytrope cannot hold 
if we remove the boundary of the adiabatic wall. As a 
demonstration, Fig|3Jc) plots the results removing the 
boundary, in which the initial state is the same distribu- 
tion as in run A. As is expected, the high-energy particles 
freely escape outwards from the central region and the re- 
sultant distribution function /(e) sharply falls off at the 
energy region £ ~ 0, indicating that the density contrast 
D becomes effectively large. Thus, compared to the sys- 
tem confined in the wall, the removal of the boundary 
makes the stellar system unstable and the core-collapse 
takes place earlier. Nevertheless, focusing on the inner 
part of the denser region, the evolution of the core is not 
significantly affected by the escape particles at the outer 
part and can be fitted by a sequence of stellar polytropes 
(see also the dashed line in Figl^^c)). The successful fit to 
the density profile p{r) almost remains the same. Hence, 
the stellar polytrope characterized by the Tsallis entropy 
can be even realized in a realistic situation removing the 
boundary of the adiabatic wall. 

Summary & Discussions. — We have performed a set 
of numerical simulation of long-term stellar dynamical 
evolution away from the BG state and found that the 
transient state of the system confined in an adiabatic 
wall can be remarkably fitted by a sequence of stellar 
polytropes. This is even true in the case removing the 
outer boundary. Therefore, the stellar polytropic distri- 



bution can be a quasi-attractor and a quasi-equilibrium 
state of a self-gravitating system. 

Alternative characterization of the transients away 
from the BG state might be possible besides the q- 
exponential distribution of stellar polytropes. For an 
empirical characterization of observed structure, the one- 
parameter family of truncated exponential distributions, 
so-called King model has been used in the literature P, 
0,123. Also, the sequence of King model has been found 
to characterize the evolutionary sequence of density pro- 
file for isolated stellar systems without boundary We 
have also tried to fit the simulation data to the King 
model. Similarly to the stellar polytrope, the King model 
accurately describes the simulated density profile p{r) 
confined in an adiabatic wall, however, it fails to match 
the simulated distribution function /(e), especially at the 
cutoff energy scales. Therefore, from the quantitative de- 
scription of the entire phase-space structure, the power- 
law distribution of the stellar polytropes can be a better 
characterization of the quasi-equilibrium state and this 
could yield an interesting explanation of the origin of the 
empirical King model. 
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